<html lang="en">
<head>
<meta http-equiv="Content-Type" content="text/html; charset=utf-8">
<title>mpt_demo1</title>
<style type="text/css">
	body {background-color: white; color: black; font-family:sans-serif; font-size:medium;}
	a:link {color: #3300ff;}
	a:visited {color: #663399;}
	a:hover {color:#0099ff;}
	a:active {color: #0066cc;}
	a.button {text-decoration:none;}
	
	table.nav  {background-color: #dbddff;}
	table.body {margin-top:2ex; margin-bottom:2ex;}
	table.programlistingindent {margin-left:32px;}
	
	img { margin-bottom:0px; margin-top:0px;}
	tt {margin-left:0.5em; margin-right:0.5em; font-weight:lighter;}
	
	p {margin-top:0ex;}
	p.synopsis {margin-left:32px;}
	p.programlistingindent {margin-left:32px;}
	p.citetitle {margin-left:2em;}
	
	ul ul {list-style-type:square;}
	ul li p {margin-top:0ex; margin-bottom:.5ex; padding:0}
	ol li p {margin-top:0ex; margin-bottom:.5ex; padding:0}
	
	h1.reftitle {color:#a90000;}
	h1.reftitle {font-size:3.7ex; margin-top:0; margin-bottom:0; font-weight:bold}
	h1.title {color:black; font-size:4ex; margin-top:1ex; font-weight:bold}
	h2.title {color:#bd0000; margin-top:1ex; margin-bottom:.9ex; font-weight:bold; font-size:3ex}
	h3.title {color:#bd0000; margin-top:1ex; margin-bottom:.9ex; font-weight:bold; font-size:2.5ex}
	h4.title {color:#bd0000; margin-top:1ex; margin-bottom:.9ex; font-weight:bold; font-size:2ex}
	h2 {color:#bd0000; margin-top:1ex; margin-bottom:.9ex; font-weight:bold; font-size:2.5ex}
	h3 {color:#bd0000; margin-top:1ex; margin-bottom:.9ex; font-weight:bold; font-size:2ex} 
	
	pre.programlisting {margin-left:32px;}
	pre.synopsis {margin-left:32px;}
	
	
	.categorytitle {margin-top:8px; padding-top:0px;}
	.categorylist {background-color: #e1e6f2;}
 	</style>
</head>
<body>
<a name="top_of_page"></a><p style="font-size:1px;"></p>
<h1 class="reftitle">mpt_demo1</h1>
<h2>Purpose</h2>
<p>Demonstration of basic usage of the geometric library</p>
<h2>Syntax</h2>
<pre class="synopsis">mpt_demo1</pre>
<h2>Description</h2>
<p></p>
        Basic usage of the new interface to geometric library
	<h2>Example(s)</h2>
<h3>Example 
				1</h3>Create polytope:<pre class="programlisting">P = Polyhedron('V', randn(10,2)); </pre>
<pre class="programlisting"></pre> Plot the polytope <pre class="programlisting"> P.plot('color','b','alpha',0.3,'linewidth',1,'linestyle','--'); axis off; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_1.png" alt="../../../fig/mpt/demos/mpt_demo1_img_1.png" width="60%"></p> Double description created and stored automatically <pre class="programlisting"> P </pre>
<pre class="programlisting">Polyhedron in R^2 with representations:
    H-rep (irredundant) : Inequalities   4 | Equalities   0
    V-rep (irredundant) : Vertices   4 | Rays   0
Functions : none
</pre>
<h3>Example 
				2</h3> Create Polyhedron <pre class="programlisting">P = Polyhedron('V', randn(10,2), 'R', randn(1,2)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting">P.plot;</pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_2.png" alt="../../../fig/mpt/demos/mpt_demo1_img_2.png" width="60%"></p> Create another Polyhedron <pre class="programlisting"> P = Polyhedron('V', randn(50,3), 'R', randn(1,3)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> P.plot;  axis vis3d; figure(1); </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_3.png" alt="../../../fig/mpt/demos/mpt_demo1_img_3.png" width="60%"></p>
<h3>Example 
				3</h3> Lower-dimensional polyhedra <pre class="programlisting"> P = Polyhedron('H', [randn(30,3) ones(30,1)], 'He', [randn(1,3) 0]); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> P.plot('alpha',0.3); </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_4.png" alt="../../../fig/mpt/demos/mpt_demo1_img_4.png" width="60%"></p>
<h3>Example 
				4</h3> Boxes  <pre class="programlisting"> P = Polyhedron('lb', -rand(3,1), 'ub', rand(3,1), 'H', [randn(10,3) 1.5*ones(10,1)]); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> P.plot('alpha',0.3); </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_5.png" alt="../../../fig/mpt/demos/mpt_demo1_img_5.png" width="60%"></p>
<h3>Example 
				5</h3> Polyhedron queries  Project onto polyhedron: <pre class="programlisting"> P = Polyhedron(randn(10,2), ones(10,1)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> P.plot('color','g'); 
                  hold on;
                  axis square;            
                  x = 5*randn(2,1);
                  sol = P.project(x); 
                  pplot([x sol.x]','bo');
                  sep = P.separate(x);
                  v = axis; 
                  s = Polyhedron('He', sep, 'lb', [v(1);v(3)],'ub',[v(2);v(4)]); 
                  s.plot; 
                  sol = P.interiorPoint; 
                  pplot(sol.x, 'ro'); 
            </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_6.png" alt="../../../fig/mpt/demos/mpt_demo1_img_6.png" width="60%"></p> Addition <pre class="programlisting"> P = Polyhedron(3*randn(10,2)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> Q = Polyhedron('H', [randn(10,2) ones(10,1)]); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> P.plot; hold on; Q.plot; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_7.png" alt="../../../fig/mpt/demos/mpt_demo1_img_7.png" width="60%"></p>
<pre class="programlisting"> PQ = P+Q; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> PQ.plot('alpha',0.1,'color','r'); </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_8.png" alt="../../../fig/mpt/demos/mpt_demo1_img_8.png" width="60%"></p>Lower-dimensional addition <pre class="programlisting">
                for i=1:5
                    P(i) = Polyhedron('V', [0 0 0;randn(1,3)]);
                end
            </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> plot(P,'linewidth',2); hold on;
                Q = Polyhedron;
                for i=1:5
                Q = Q + P(i);
                end
                Q.plot('alpha',0.1,'color','b');
                axis vis3d
            </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_9.png" alt="../../../fig/mpt/demos/mpt_demo1_img_9.png" width="60%"></p>
<h3>Example 
				6</h3> Operations with convex sets <pre class="programlisting"> x = sdpvar(2,1); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> T = eye(2); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> F = [x &lt;= 0.1 ; x'*T'*T*x &lt;= 1]; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> Y = YSet(x, F); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> Y.plot('alpha',0.3,'color','b','grid',50); 
                hold on;
                x = [-2;-2];
                sol = Y.project(x);
                pplot([x sol.x]','bo');
                sol = Y.separate(x); 
                v = axis;
                s = Polyhedron('He', sol, 'lb', [v(1);v(3)],'ub',[v(2);v(4)]);
                s.plot('alpha',0.2);
                O = Y.outerApprox;
                O.plot('alpha',0.1,'color','b'); axis(axis*1.1);
            </pre>
<pre class="programlisting">Plotting...
35 of 50
</pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_10.png" alt="../../../fig/mpt/demos/mpt_demo1_img_10.png" width="60%"></p>
<h3>Example 
				7</h3>Create set of polytopes<pre class="programlisting">Ps = PolyUnion; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting">
                for i=1:5
                  Ps.add(Polyhedron(randn(10,2)) + 5*randn(2,1));
                end
            </pre>
<pre class="programlisting"></pre>
<pre class="programlisting">
                for i=1:5
                   Ps.add(Polyhedron('H', [randn(10,2) ones(10,1)], 'He', [randn(1,2) 0]) + 5*randn(2,1));
                end
            </pre>
<pre class="programlisting"></pre>
<pre class="programlisting">
                for i=1:5
                   Z = Polyhedron;
                   for j=1:3,
                      Z = Z + Polyhedron([0 0;randn(1,2)]);
                   end
                   Ps.add(Z+5*randn(2,1));
                end
            </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> Ps.plot; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_11.png" alt="../../../fig/mpt/demos/mpt_demo1_img_11.png" width="60%"></p> Create complex <pre class="programlisting"> T = Polyhedron('H', [randn(30,2) ones(30,1)]); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> T = triangulate(T); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> T.plot; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_12.png" alt="../../../fig/mpt/demos/mpt_demo1_img_12.png" width="60%"></p>
<h3>Example 
				8</h3> Polyhedral functions <pre class="programlisting"> F = Polyhedron('V', 2*randn(10,2))</pre>
<pre class="programlisting">Polyhedron in R^2 with representations:
    H-rep               : Unknown (call computeHRep() to compute)
    V-rep (redundant)   : Vertices  10 | Rays   0
Functions : none
</pre>
<pre class="programlisting"> F.addFunction(Function(@(x) sin(x(1))*cos(x(2))),'func1'); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> F.fplot; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_13.png" alt="../../../fig/mpt/demos/mpt_demo1_img_13.png" width="60%"></p> PolyUnion <pre class="programlisting"> func = @(x) sin(x(1))*cos(x(2)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> 
                for i=1:5
                  F(i) = Polyhedron('V', 2*randn(10,2)) + 5*randn(2,1);
                end
                F.addFunction(Function(@(x) sin(x(1))*cos(x(2))),'wave');
            </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> Ps = PolyUnion(F); </pre>
<pre class="programlisting"></pre> Plot the function over the set <pre class="programlisting"> Ps.fplot; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_14.png" alt="../../../fig/mpt/demos/mpt_demo1_img_14.png" width="60%"></p>
<h3>Example 
				9</h3> Parametric solutions  Formulate MPC problem with the following dimensions <pre class="programlisting"> n = 2; m = 1; N = 5; </pre>
<pre class="programlisting"></pre> Setup the process model and constraints <pre class="programlisting"> A = [1 1; 0 1]; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> B = [1; 0.5]; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> xlb = -10*ones(n,1); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> xub = 10*ones(n,1); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> ulb = -5*ones(m,1); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> uub = 5*ones(m,1); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> R = 2*eye(m); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> Q = 0.2*eye(n); </pre>
<pre class="programlisting"></pre>Formulate MPQP using Yalmip <pre class="programlisting"> x = sdpvar(n,N,'full'); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> u = sdpvar(m,N-1,'full'); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> cost = 0; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> F = []; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> 
                    for i=1:N-1
                        F = F + (x(:,i+1) == A*x(:,i) + B*u(:,i));
                        F = F + (xlb &lt;= x(:,i) &lt;= xub);
                        F = F + (ulb &lt;= u(:,i) &lt;= uub);  
                        if i &gt; 1
                            cost = cost + x(:,i)'*Q*x(:,i);
                        end
                        cost = cost + u(:,i)'*R*u(:,i);
                    end
                </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> F = F + [xlb &lt;= x(:,end) &lt;= xub]; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> cost = cost + x(:,end)'*Q*x(:,end);</pre>
<pre class="programlisting"></pre> Solve using MPQP solver  Change globally the parametric QP solver <pre class="programlisting"> mptopt('pqpsolver','MPQP'); </pre>
<pre class="programlisting"></pre> Construct the problem <pre class="programlisting"> problem1 = Opt(F, cost, x(:,1), u(:)); </pre>
<pre class="programlisting"></pre> Solve <pre class="programlisting"> res1 = problem1.solve; </pre>
<pre class="programlisting">Calling mpt_mpqp_26 with default options...
mpt_mpqp: 13 regions
</pre> Solve using PLCP solver  Change globally the parametric QP solver <pre class="programlisting"> mptopt('pqpsolver','PLCP'); </pre>
<pre class="programlisting"></pre> Call problem constructor <pre class="programlisting"> problem2 = Opt(F, cost, x(:,1), u(:)); </pre>
<pre class="programlisting"></pre> Solve <pre class="programlisting"> res2 = problem2.solve; </pre>
<pre class="programlisting">mpt_plcp: 13 regions
</pre> Plot the partitions <pre class="programlisting"> plot(res1.xopt);  title('PLCP'); axis tight; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_15.png" alt="../../../fig/mpt/demos/mpt_demo1_img_15.png" width="60%"></p>
<pre class="programlisting"> plot(res2.xopt);  title('MPQP'); axis tight; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_16.png" alt="../../../fig/mpt/demos/mpt_demo1_img_16.png" width="60%"></p> Formulate MPLP problem in YALMIP<pre class="programlisting">cost = 0; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> F = []; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting">
                for i=1:N-1
                    F = F + (x(:,i+1) == A*x(:,i) + B*u(:,i));
                    F = F + (xlb &lt;= x(:,i) &lt;= xub);
                    F = F + (ulb &lt;= u(:,i) &lt;= uub);    
                    if i &gt; 1
                        cost = cost + norm(Q*x(:,i),1);
                    end
                    cost = cost + norm(R*u(:,i),1);
                end
            </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> F = F + [xlb &lt;= x(:,end) &lt;= xub]; </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> cost = cost + norm(Q*x(:,end),1); </pre>
<pre class="programlisting"></pre> Solve using MPLP solver <pre class="programlisting"> mptopt('plpsolver','MPLP'); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> problem3 = Opt(F, cost, x(:,1), u(:)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> res3 = problem3.solve; </pre>
<pre class="programlisting">Calling mpt_mplp_26 with default options...
mpt_mplp: 28 regions
</pre> Solve using PLCP solver <pre class="programlisting"> mptopt('plpsolver','PLCP'); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> problem4 = Opt(F, cost, x(:,1), u(:)); </pre>
<pre class="programlisting"></pre>
<pre class="programlisting"> res4 = problem4.solve; </pre>
<pre class="programlisting">regions:   21, unexplored: 5 
mpt_plcp: 28 regions
</pre> Plot the partitions <pre class="programlisting"> plot(res3.xopt);  title('MPLP'); axis tight; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_17.png" alt="../../../fig/mpt/demos/mpt_demo1_img_17.png" width="60%"></p>
<pre class="programlisting"> plot(res4.xopt);  title('PLCP'); axis tight; </pre>
<pre class="programlisting"></pre>
<p class="programlistingindent"><img src="../../../fig/mpt/demos/mpt_demo1_img_18.png" alt="../../../fig/mpt/demos/mpt_demo1_img_18.png" width="60%"></p>
<h2>See Also</h2>
<a href="./mpt_demo_sets1.html">mpt_demo_sets1</a>, <a href="./mpt_demo_functions1.html">mpt_demo_functions1</a>, <a href="./mpt_demo_unions1.html">mpt_demo_unions1</a><p></p>
<table class="nav" summary="Navigation aid" border="0" width="100%" cellpadding="0" cellspacing="0"><tr valign="top">
<td align="left" width="20">
<a href="mpt_demo2.html" class="button">&#9664;</a>  </td>
<td align="left">mpt_demo2</td>
<td>  </td>
<td align="right">mpt_demo_unions2</td>
<td align="right" width="20"><a href="mpt_demo_unions2.html" class="button">&#9654;</a></td>
</tr></table>
<br><p>©  <b>2010-2013</b>     Colin Neil Jones: EPF Lausanne,    <a href="mailto:colin.jones@epfl.ch">colin.jones@epfl.ch</a></p>
</body>
</html>
